# This script reproduces Figure 10 #

rm(list=ls())
options(digits=4)
pkg <- c("PanelMatch", "ggplot2", "plm", 
         "matrixStats",
         "foreign", 
         "stringi", "stringdist",
         "DataCombine","lmtest", "multiwayvcov",
         "data.table",
         "dplyr", 
         "stargazer",
         "clubSandwich", "mediation",
         "doBy")
lapply(pkg, require, character.only = TRUE)

# load datasets
proprietary <- read.csv("proprietary.csv")
load("non_p.RData")
source("functions.R")
# merge 
d_sub1 <- merge(d_sub1, proprietary, by = c("county_id", "time"), all.x = T)


### county-month, land
set.seed(2021)
times <- unique(d_sub1$time)
reps <- 500
land_points <- rep(NA, reps)
land_tstats <- rep(NA, reps)
for (j in 1:reps){
  d_sub1$tour_alt <- NA
  for (i in 1:length(times)) {
    d_sub1$tour_alt[d_sub1$time == times[i]] <- sample(x = d_sub1$tour[d_sub1$time == times[i]], 
                                                       size = length(d_sub1$tour[d_sub1$time == times[i]]))
  }

  formula <- as.formula(paste0(paste0("log_transaction_area", "_f", 1), 
                               "~ tour_alt + 
                               log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                               log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                               log_expd_l1 + log_middle_school_l1 + 
                               # msec2currentsec2_l1 +      
                               msec2currentsec_l1 + 
                               msec2currentgvn_l1 + 
                               mayor2currentsec_l1 + 
                               #mayor2currentsec2_l1 + 
                               mayor2currentgvn_l1 + 
                               as.factor(county_id) + as.factor(time)"))
  
  model_land_con <- lm(formula = formula,

                       data = d_sub1)
  vcov_land_con <- clubSandwich::vcovCR(model_land_con, type="CR1S", cluster = d_sub1$county_id)
  land_points[j] <- coeftest(model_land_con, vcov_land_con)["tour_alt", 1]
  land_tstats[j] <- coeftest(model_land_con, vcov_land_con)["tour_alt", 3]
  cat("The", j, "th value is",  land_points[j], "\n")
}


formula <- as.formula(paste0(paste0("log_transaction_area", "_f", 1), 
                             "~ tour + 
                             log_transaction_area_l1 + log_transaction_area_l2 + log_transaction_area_l3 + 
                             log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                             log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                             log_expd_l1 + log_middle_school_l1 + 
                             # msec2currentsec2_l1 +      
                             msec2currentsec_l1 + 
                             msec2currentgvn_l1 + 
                             mayor2currentsec_l1 + 
                             #mayor2currentsec2_l1 + 
                             mayor2currentgvn_l1 +
                             as.factor(county_id) + as.factor(time) "))
model_land_con <- lm(formula = formula,
                     data = d_sub1)
vcov_land_con <- clubSandwich::vcovCR(model_land_con, type="CR1S", cluster = d_sub1$county_id)
land_point <- coeftest(model_land_con, vcov_land_con)["tour", 1]
land_tstat <- coeftest(model_land_con, vcov_land_con)["tour", 3]

land_placebo_results <- list("land_points" = land_points, 
                             "land_tstats" = land_tstats,
                             "land_point" = land_point,
                             "land_tstat" = land_tstat)
save(land_placebo_results, 
     file = "output/land_placebo_results.RData")
load("output/land_placebo_results.RData")
pdf(file = "output/Figure10.pdf", 
    width = 10, height = 8)
par(mfrow = c(1,2))
hist(land_placebo_results$land_tstats,
     main = "",
     breaks = 40,
     xlab = "T-statistics for Land Auction Results"
)
abline(v = land_placebo_results$land_tstat, col = "red",
       lty = "dashed")
text(land_placebo_results$land_tstat + 0.8, 40, 
     paste0("p = ", 
            sum(land_placebo_results$land_tstats <  land_placebo_results$land_tstat)/500
     ), col = "blue")

hist(land_placebo_results$land_points,
     main = "",
     breaks = 40,
     xlab = "Point Estimates for Land Auction Results"
)
abline(v = land_placebo_results$land_point, col = "red",
       lty = "dashed")
text(land_placebo_results$land_point+.05, 31.8, 
     paste0("p = ", 
            sum(land_placebo_results$land_points <  land_placebo_results$land_point)/500
     ), col = "blue")
dev.off()
